Modeling Corn Price Series
Corn Futures Price
In the following, we analyze the Continuous Futures
Contract for Corn between
01/01/2003 to 12/31/2020.
Find more about Quandl
Continuous Futures and Corn
Futures
The price is in USD/bushel.
In the plot of the log-returns we can clearly see a lot of volatility clusters of different size and persistence over the time period we are analyzing; this is confirmed in the ACF plot of the squared series of the log-returns which shows that they are autoccorrelated.
Normal ARMA-GARCH
To fit an \(ARMA(p,q)-GARCH(p_g,q_g)\) model we have a standard ARMA(p,q)-model of the form \[\phi(B)(r_t-\mu)=\theta(B)X_t\] where the innovations \(X_t\) are a \(GARCH(p_g,q_g)\) process s.t. we have:
\[ \phi(B)(r_t-\mu)=\theta(B)X_t, \quad X_t=\sigma_t\epsilon_t, \quad \sigma_t^2=\alpha_0+\sum_{h=1}^{p_g}\alpha_hX_{t-h}^2+\sum_{j=1}^{q_g}\gamma_j\sigma_{t-j}^2 \]
Particularly if we fit an ARMA(1,1)-GARCH(1,1) model we have:
\[ \begin{aligned} r_t &= \mu_t + X_t ,\quad X_t = \sigma_t\epsilon_t, \quad\epsilon_t\sim N(0,1)\\ \mu_t &= \mu(1-\phi_1) + \phi_1r_{t-1} + \theta_1X_{t-1} \\ \sigma_t^2 &= \alpha_0 + \alpha_1X_{t-1}^2 + \gamma_1\sigma_{t-1}^2 \end{aligned} \]
Model Specification
# ------------------------------------------------------------------------------
# --- Model fit using the "rugarch" package ---
# ------------------------------------------------------------------------------
# Create the model, specify the variance.model, the mean.model and the distribution.model
# We set the orders of our model as GARCH(1,1) and ARMA(1,1)
model <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(1, 1), include.mean = TRUE),
distribution.model = "norm")
model@model$pars## Level Fixed Include Estimate LB UB
## mu 0 0 1 1 NA NA
## ar1 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## beta1 0 0 1 1 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 0 0 NA NA
## shape 0 0 0 0 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NA# Estimation/fitting of the model
model_fit <- ugarchfit(spec = model, data = lr)
# Fitted Model
model_fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : norm
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## mu -0.000229 0.000205 -1.1198 0.262809
## ar1 -0.894389 0.040326 -22.1792 0.000000
## ma1 0.911589 0.037742 24.1531 0.000000
## omega 0.000002 0.000004 0.5632 0.573300
## alpha1 0.051315 0.025554 2.0081 0.044635
## beta1 0.942541 0.028210 33.4122 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## mu -0.000229 0.000780 -0.293716 0.768975
## ar1 -0.894389 0.232442 -3.847792 0.000119
## ma1 0.911589 0.215950 4.221296 0.000024
## omega 0.000002 0.000066 0.032611 0.973985
## alpha1 0.051315 0.449528 0.114153 0.909116
## beta1 0.942541 0.493670 1.909253 0.056229
##
## LogLikelihood : 12464.75
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.5006
## Bayes -5.4921
## Shibata -5.5006
## Hannan-Quinn -5.4976
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.01409 0.9055
## Lag[2*(p+q)+(p+q)-1][5] 1.19434 0.9999
## Lag[4*(p+q)+(p+q)-1][9] 3.25805 0.8477
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.8551 0.3551
## Lag[2*(p+q)+(p+q)-1][5] 1.4526 0.7516
## Lag[4*(p+q)+(p+q)-1][9] 4.2091 0.5538
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.3945 0.500 2.000 0.5299
## ARCH Lag[5] 1.1004 1.440 1.667 0.7034
## ARCH Lag[7] 1.9341 2.315 1.543 0.7313
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 75.5507
## Individual Statistics:
## mu 0.07968
## ar1 0.11323
## ma1 0.09977
## omega 11.96273
## alpha1 0.23464
## beta1 0.24609
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.1121 0.9107
## Negative Sign Bias 0.6952 0.4870
## Positive Sign Bias 0.7350 0.4624
## Joint Effect 1.0490 0.7894
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 120.9 7.597e-17
## 2 30 153.0 8.261e-19
## 3 40 159.2 1.874e-16
## 4 50 180.0 7.424e-17
##
##
## Elapsed time : 0.492682# empirical density plot
plot(model_fit, which=8)
# normal QQ-plot
plot(model_fit, which=9)
# estimated coefficients
coef(model_fit)## mu ar1 ma1 omega alpha1
## -2.292451e-04 -8.943893e-01 9.115892e-01 2.159513e-06 5.131501e-02
## beta1
## 9.425407e-01Note:
Clearly the statistics of our fitted model suggest that many of our parameters are not statistically significant and maybe the model order is not adequate. Also notice that the Pearson gof test rejects \(H_0\) and looking at the empirical density and QQ plots we can see that the assumption of normality (\(\epsilon_t\sim N(0,1)\)) is underestimating the tails and may not be a good fit. Nevertheless, for comparison reasons we will show the ARMA-GARCH with N(0,1) fit in the following plots, and after that we will proceed to fit other model orders with with a Student-t distribution.
The parameters of the fitted ARMA(1,1)-GARCH(1,1) correspond to:
- mu (\(\mu\)): estimated mean of the ARMA process.
- ar1 (\(\phi_1\)): autoregressive parameter for the process \(r_t\)
- ma1 (\(\theta_1\)): moving average parameter for the innovations \(X_t\)
- omega (\(\alpha_0\)): estimated (constant) drift of our conditional variance model \(\sigma_t^2\).
- alpha1 (\(\alpha_1\)): estimated weight of how much our process \(\sigma_t\) depends on past values of the squared process \(X_{t-1}^2\).
- beta1 (\(\gamma_1\)): estimated weight of how much our process \(\sigma_t\) depends on past volatilities \(\sigma_{t-1}^2\)
It is worth noting that beta1: \(|\gamma_1|=0.942541<1\) and alpha1+beta1: \(|\alpha_1 + \gamma_1|=0.993869<1\) which tells us that our model is both invertible and stationary, and since the sum \(|\alpha_1+\gamma_1|\) is close to \(1\), we can see that volatility persists for a long time which is consistent with what we saw in the plot of the log-returns.
Model diagnostics (standard residuals)
The model diagnostics is done based on the standardized log-returns
\[ \begin{aligned} \hat{\epsilon_t}=\frac{r_t - \mu_t(\hat{\Theta})}{\sigma_t(\hat{\Theta})} \end{aligned} \]
Standardized residuals
- We can see that volatility clusters are no longer present in the standardized log-returns series. Nevertheless we can still see very large and small values in the log-returns so we will try specifying a Student-t distribution for the model to account for heavier tails in the log-returns.
- The ACF of \(\hat{\epsilon}_t\) shows no autocorrelation suggesting that the model for the conditional mean is appropriate.
- The ACF of \(\hat{\epsilon}_t^2\) shows almost no autocorrelation (with the exception of lags 9 and 30) suggesting that we got rid of the volatility clusters, but remember that we found out that the parameters were not significant for the volatility model so we still need to address this.
Still the main drawback is the apparent non-normality of the log-returns and we are assuming that positive and negative returns have the same effect because in the GARCH model volatility depends squared returns. To change this we could consider using an skewed version of our distribution for the model. i.e. “snorm”, “sstd”.
Student-t ARMA-GARCH
From the ARMA-GARCH model that we previously fitted, we saw that the model specification was not altogether correct (based on the t statistic).
After testing various model orders, we find that the best fit does correspond to an ARMA-GARCH model but with a Student-t model distribution.
So we will now fit an ARMA(1,1)-GARCH(1,1) (without mean \(\mu\) since we saw that it was not significant so we take it as zero) and with a Student-t distribution for the model: \(\epsilon_t\sim t_{\nu}\)
Model Specification
# Create the model, specify the variance.model, mean.model and the distribution.model
#model2 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
# mean.model = list(armaOrder = c(1, 1), include.mean = FALSE),
# distribution.model = "std", fixed.pars = list(omega=0))
model2 <- ugarchspec(variance.model = list(model = "sGARCH", garchOrder = c(1, 1)),
mean.model = list(armaOrder = c(1, 1), include.mean = FALSE),
distribution.model = "std")
mydata <- data.frame(lr)
row.names(mydata) <- row.names(corn_ts)[-c(1:2)]
model2@model$pars## Level Fixed Include Estimate LB UB
## mu 0 0 0 0 NA NA
## ar1 0 0 1 1 NA NA
## ma1 0 0 1 1 NA NA
## arfima 0 0 0 0 NA NA
## archm 0 0 0 0 NA NA
## mxreg 0 0 0 0 NA NA
## omega 0 0 1 1 NA NA
## alpha1 0 0 1 1 NA NA
## beta1 0 0 1 1 NA NA
## gamma 0 0 0 0 NA NA
## eta1 0 0 0 0 NA NA
## eta2 0 0 0 0 NA NA
## delta 0 0 0 0 NA NA
## lambda 0 0 0 0 NA NA
## vxreg 0 0 0 0 NA NA
## skew 0 0 0 0 NA NA
## shape 0 0 1 1 NA NA
## ghlambda 0 0 0 0 NA NA
## xi 0 0 0 0 NA NA# Estimation/fitting of the model
model2_fit <- ugarchfit(spec = model2, data = mydata, solver = "solnp")
# Fitted Model
model2_fit##
## *---------------------------------*
## * GARCH Model Fit *
## *---------------------------------*
##
## Conditional Variance Dynamics
## -----------------------------------
## GARCH Model : sGARCH(1,1)
## Mean Model : ARFIMA(1,0,1)
## Distribution : std
##
## Optimal Parameters
## ------------------------------------
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.908924 0.010598 -85.7610 0.000000
## ma1 0.925411 0.011055 83.7119 0.000000
## omega 0.000002 0.000001 1.9910 0.046484
## alpha1 0.055366 0.007431 7.4505 0.000000
## beta1 0.939458 0.008057 116.6040 0.000000
## shape 5.299195 0.374311 14.1572 0.000000
##
## Robust Standard Errors:
## Estimate Std. Error t value Pr(>|t|)
## ar1 -0.908924 0.010186 -89.23405 0.000000
## ma1 0.925411 0.007651 120.94891 0.000000
## omega 0.000002 0.000003 0.72166 0.470504
## alpha1 0.055366 0.026507 2.08873 0.036732
## beta1 0.939458 0.027128 34.63115 0.000000
## shape 5.299195 0.666284 7.95336 0.000000
##
## LogLikelihood : 12654.72
##
## Information Criteria
## ------------------------------------
##
## Akaike -5.5844
## Bayes -5.5759
## Shibata -5.5844
## Hannan-Quinn -5.5814
##
## Weighted Ljung-Box Test on Standardized Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.4022 0.5260
## Lag[2*(p+q)+(p+q)-1][5] 0.7405 1.0000
## Lag[4*(p+q)+(p+q)-1][9] 2.1688 0.9762
## d.o.f=2
## H0 : No serial correlation
##
## Weighted Ljung-Box Test on Standardized Squared Residuals
## ------------------------------------
## statistic p-value
## Lag[1] 0.04003 0.8414
## Lag[2*(p+q)+(p+q)-1][5] 0.26910 0.9868
## Lag[4*(p+q)+(p+q)-1][9] 0.65095 0.9964
## d.o.f=2
##
## Weighted ARCH LM Tests
## ------------------------------------
## Statistic Shape Scale P-Value
## ARCH Lag[3] 0.1982 0.500 2.000 0.6562
## ARCH Lag[5] 0.2248 1.440 1.667 0.9595
## ARCH Lag[7] 0.4918 2.315 1.543 0.9791
##
## Nyblom stability test
## ------------------------------------
## Joint Statistic: 74.9014
## Individual Statistics:
## ar1 0.10170
## ma1 0.08659
## omega 8.88725
## alpha1 0.77875
## beta1 0.65982
## shape 0.61138
##
## Asymptotic Critical Values (10% 5% 1%)
## Joint Statistic: 1.49 1.68 2.12
## Individual Statistic: 0.35 0.47 0.75
##
## Sign Bias Test
## ------------------------------------
## t-value prob sig
## Sign Bias 0.1189 0.9054
## Negative Sign Bias 0.3799 0.7040
## Positive Sign Bias 0.3919 0.6952
## Joint Effect 0.3274 0.9548
##
##
## Adjusted Pearson Goodness-of-Fit Test:
## ------------------------------------
## group statistic p-value(g-1)
## 1 20 12.03 0.8843
## 2 30 17.15 0.9598
## 3 40 33.69 0.7103
## 4 50 28.12 0.9928
##
##
## Elapsed time : 0.7097969# empirical density
plot(model2_fit, which=8)
# t-student QQ-plot
plot(model2_fit, which=9)
# estimated coefficients
coef(model2_fit)## ar1 ma1 omega alpha1 beta1
## -9.089240e-01 9.254107e-01 2.090566e-06 5.536554e-02 9.394575e-01
## shape
## 5.299195e+00The parameters of the fitted Student-t ARMA(1,1)-GARCH(1,1) correspond to:
- ar1 (\(\phi_1\)): autoregressive parameter for the process \(r_t\)
- ma1 (\(\theta_1\)): moving average parameter for the innovations \(X_t\)
- omega (\(\alpha_0\)): estimated (constant) drift of our conditional variance model \(\sigma_t^2\).
- alpha1 (\(\alpha_1\)): estimated weight of how much our process \(\sigma_t\) depends on past values of the squared process \(X_{t-1}^2\).
- beta1 (\(\gamma_1\)): estimated weight of how much our process \(\sigma_t\) depends on past volatilities \(\sigma_{t-1}^2\)
- shape (\(\nu\)): estimated degrees of freedom for the Student-t distribution.
Again we have beta1: \(|\gamma_1|=0.9394575<1\) and alpha1+beta1: \(|\alpha_1 + \gamma_1|=0.999<1\) so our model is both invertible and stationary and the volatility persists for a long time which is consistent with our log-returns \(r_t\).
This time notice how the Student-t distribution fits very well out data as shown in both the empirical density and QQ plots.
We can see that our confidence bands (at the 95% C.I.) cover the log-returns very well in line with the volatility clusters of the series; the exceedances are approx <5% according to our C.I. and only in a few cases do they exceed by a considerable amount out bands.
Model diagnostics (standard residuals)
Standardized residuals
- We can see that volatility clusters are no longer present in the standardized log-returns series.
- The ACF of \(\hat{\epsilon}_t\) shows no autocorrelation suggesting that the model for the conditional mean is appropriate.
- The ACF of \(\hat{\epsilon}_t^2\) shows no autocorrelation suggesting that we got rid of the volatility clusters and that the order of the model is appropriate.